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The effect of surface disorder on electronic systems is particularly interesting for topological 
phases with surface and edge states. Using exact diagonalization, it has been demonstrated that the 
surface states of a 3D topological insulator survive strong surface disorder, and simply get pushed 
to a clean part of the bulk. Here we explore a new method which analytically eliminates the clean 
bulk, and reduces a D-dimensional problem to a Hamiltonian-diagonalization problem within the 
{D — /-dimensional disordered surface. This dramatic reduction in complexity allows the analysis 
of significantly bigger systems than is possible with exact diagonalization. We use our method to 
analyze a 2D topological spin-Hall insulator with non-magnetic and magnetic edge impurities, and 
we calculate the probability density (or local density of states) of the zero-energy eigenstates as 
a function of edge-parallel momentum and layer index. Our analysis reveals that the system size 
needed to reach behavior in the thermodynamic limit increases with disorder. We also compute the 
edge conductance as a function of disorder strength, and chart a lower bound for the length scale 
marking the crossover to the thermodynamic limit. 


I. INTRODUCTION 

Solid state systems inevitably contain impurities. 
Study of the impurity effects is of special importance in 
topological insulators (for a review see, e.g. Ref. [THU), be¬ 
cause their surface states are expected to be fundamen¬ 
tally robust with respect to certain types of disorder. Un¬ 
derstanding how impurities affect these systems also has 
practical implications because the fabrication of topolog¬ 
ical insulators (TIs) requires fine tuning of the doping 
concentration to place the Fermi energy within the band 
gajP^^ or to induce superconductivitj^^®^. Impurities 
that are located at or near the TI surface are especially 
interesting. One relevant example is the treatment with 
NO2 (see supplement of Ref. lU on the surface of as-grown 
topological insulator Bi2_a;Cua;Se3 which is necessary to 
prevent the surface band bending caused by the adsorp¬ 
tion of residual atoms present in the vacuum chamber. 
As the manifestations of topological systems - such as 
the protected surface states or Majorana fermions - are 
localized near the surface or the edge of the systenJiTE^, 
the effect of surface impurities on the surface spectral and 
transport properties are of both theoretical and practical 
interespSHi^. Local density of states (LDOS) can be effi¬ 
ciently probed by scanning tunneling spectroscopy while 
the surface spectral function can be extracted from angle 
resolved photoemission spectroscopy (ARPES) measure¬ 
ments. Various techniques have been employed to study 
the transport properties of TI surfaces. 

Recently, an exact diagonalization (ED) approach 
has be en ap plied to non-interacting TIs with surface 
disordeiffSHSI. This study found evidence for a crossover 
between a nearly ballistic response of the surface elec¬ 
trons at weak disorder, localization physics at intermedi¬ 


ate disorder, and then a restored nearly ballistic surface 
state hiding in the second layer when disorder is very 
strong. Further investigations of this phenomenon using 
ED are going to prove very challenging: the computa¬ 
tional cost for analyzing the surface physics can be very 
high because the bulk of the system—although gapped— 
must be treated on equal footing with the surface de¬ 
grees of freedom. The system sizes accessible to ED 
analysis are, therefore, rather small. In particular, we 
demonstrate here that when strong surface impurities are 
present in the system, the lower bound of system size re¬ 
quired to clearly resolve the bulk electronic properties 
and their effect on the surface states becomes large in 
proportion to the impurity potential strength (at least in 
two-dimensions). Treating sufficiently large systems us¬ 
ing ED becomes computationally challenging in this limit 
and new techniques to address the problem are required. 
In this manuscript we develop such a technique. We also 
emphasize the need to treat sufficiently large systems in 
order to distinguish system properties at the thermody¬ 
namic limit versus those of the “quantum dot” regime, 
where finite size effects dominate. 

In this manuscript, we introduce a new technique that 
allows us to efficiently extract the surface state properties 
of surface-disordered TIs. We obtain properties such as 
the surface spectral function, LDOS and transport prop¬ 
erties of the surface channels, by essentially “integrating 
out” the clean bulk degrees of freedom analytically and 
obtain the effective surface-Hamiltonian describing a TI 
surface with arbitrarily strong impurities. This approach 
not only allows us to reduce the computational difficulty 
by one dimension (e.g. for a 3 D TI with surface disorder 
we only need to solve a 2 D problem), but also allow us to 
map a strong disordered problem into a weak disordered 
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one where perturbation theory is valid. By constructing 
a self-consistent transfer-matrix approach, we are able 
to recover the exact energies and wavefunctions of the 
surface states both at the disordered layer and in the 
remaining bulk layers. 


The manuscript is organized as follows. In Sec. |ITj 
we first explain how to integrate out the clean bulk de¬ 
grees of freedom. In Sec. |III| we then introduce a generic 
model Hamiltonian of a 2D TI with impurities on one of 
its edges. While our method is applicable to any layered 
system in arbitrary dimension, we choose to concentrate 
on the 2D case since it is both simple to present and 
analyze. Next, we report a series of results that clearly 
differentiate the OD quantum-dot regime from the bulk 
regime (where finite-size errors are suppressed) and spec¬ 
ify the lower bound on the system size to observe the 
latter for magnetic and non-magnetic edge impurities. 
We conclude the section with a discussion of the surface 
properties. Lastly, in Sec. IV the conductance through 
the 2D TI edge channels is computed. The latter pro¬ 
vides complementary information to the spectral proper¬ 
ties discussed in Sec. El We conclude with a discussion 
in Sec. 0 



FIG. 1. A pictorial description of a system with impurities 
on the surface only. We are interested in studying the local 
density of surface state as the impurity strength increases. 
Our strategy is to decompose the system into clean layers, 
coupled through matrices B and , and a surface containing 
impurities, and then analytically integrate out the former. 


II. GENERAL FRAMEWORK: EFFECTIVE 
SINGLE LAYER HAMILTONIAN 


A. Layered Schrodinger equation and self energy 


We begin our analysis with the Schrodinger equation 
for layers parallel to the disordered surface 

Bipn-i + [Hq + Mmp Sn,l - Ejlpn + = 0, (1) 

where ipo = 0, ipn is a wavefunction on nth layer parallel 
to the impurity surface, is an in-layer Hamiltonian, 
Vinip is an impurity potential in the first layer, B and 
Hi are hopping terms between layers, n > 0 is the layer 
index. For notational convenience, we set '0o = 0. For 
n = N, the last layer of the system, we can write exactly 

BtpN-i + [Ho - = 0 ( 2 ) 


Using the Schrodinger equation for the n = N — 1 layer, 
BlpN-2 + [Hq — Ejljjpf-i + B^-ipN = 0, and substituting 
Eq. §, we can “integrate” out the last (V’th) layer 


B'iIjn-2 + 


Ho-E + B'' 


E-Ho 


-B 


V'AT-i = 0. (3) 


Eliminating V'w introduces for ipN-i the effective poten¬ 
tial Eat-i = B^ B. By repeating this process, we 
can integrate out all layers up to the first layer and the 
following recursion relation can be found, 




E - Ho 


-B 


n+1 


( 4 ) 


with a boundary condition E^v = 0. Recall that B^ is 
a hopping to the next layer and i? is a hopping to the 
prior layer. And the effective potential E„ is obtained 
by sandwiching the Green’s function in (n -|- l)th layer 
by B'^ and B, describing a scattering process of hopping 
to the next layer, propagating, and hopping back to the 
original layer. 

Let us next write an effective Hamiltonian in the top 
layer in the following way: 


[E-Ho- EilV'i = fyn.pV’i (5) 

=U,„pV'i,, (6) 

where the recursion relation (|^ is again used to further 
simplify the clean part of the Hamiltonian. In the next 
section, we introduce a way to solve the recursion relation 
exactly. 


In this section, we provide a general derivation of our 
approach. Our goal is to exactly reduce diagonalization 
of a D-dimensional system with surface impurities to 
the diagonalization of an effective Hamiltonian describ¬ 
ing just the D — 1 dimensional surface. Taking the top 
surface to be disordered, we introduce a way to integrate 
out the clean layers from the bottom all the way to the 
top layer. This leaves us with a single layer effective 
Hamiltonian which includes the impurity potential and a 
self energy which accounts for the entire clean bulk. 


B. ‘Holographic’ mapping of the self energy 

The recursion relation can be straightforwardly solved 
by mapping the effective potential to a matrix M which 
obeys the same Schrodinger equation as the layer wave- 
functions 

+ [Ho - E]Mn + B^Mn+i = 0, (7) 

where the matrix Afy has the same dimension as the 
Hamiltonian 77||, and it is invertible by construction. 
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With Eq. Q, the recursion relation for the self-energy 
is easily solved 

= ( 8 ) 

One can directly verify that this is a solution of the re¬ 
cursion relation for satisfying boundary condition 
Mm+i = 0. For a clean bottom-surface, we can exactly 
construct Mn for a system with finite thickness (the cal¬ 
culation is detailed in the Appendix). 

The last step involves writing a close-form equation for 
the wavefunction of the top (disordered) layered. Note 
that as Mn is also a solution of the Schrodinger equation, 
an element of M„ scales with eigenvalues of transfer ma¬ 
trix of Schrodinger equation (l7|: ~ pf ”. Using 

the exact expression of M^s, we construct the left side 
of Eq. Q. Then, we obtain the Schrodinger equation 
expressed in terms of M„’s 

[BMoMf 1] V-i = MmpV'i- (9) 

This is the effective single layer Hamiltonian. The left- 
hand side contains only elements from the clean part of 
the Hamiltonian, and involves the self energy from all 
subsequent layers; the right-hand side is simply the sur¬ 
face impurity potential operating on the top-layer wave- 
function. Mq and Ml are a function of energy, and one 
can find all eigenvalues of a system by finding the ener¬ 
gies that satisfy det [BMqM^^ — Ump] = 0. The surface 
wavefunction can be subsequently found from Eq. 
which is identical to the result from exact diagonaliza- 
tion. To obtain whole wavefunction in the layers beneath 
the top layer, we apply the transfer matrix which is also 
obtained in terms of M„ as shown in the next section. 


C. Transfer matrix of wavefunctions 

The first layer wavefunction can be exactly obtained 
from Eq. therefore, the computational complexity 
is essentially reduced by one dimension. To obtain a 
full profile of the wavefunction in the subsequent layers, 
we construct an approach similar to the transfer matrix 
approach in this section. We will use the term “transfer 
matrix” quite liberally in what follows. 

E„ plays the role of effective potential in the nth layer, 
induced by integrating out the {n + l)th layer up to Nth 
layer. In the Schrodinger equation Q, such a contri¬ 
bution is accounted by the third term on the left side. 
Therefore, we have the following equality: 

B^t/Jn+I = EnV'n- (10) 

One can explicitly show this relation by the elimination 
method introduced in Eq. ([^. The transfer matrix is 
conveniently expressed in terms of the M„’s using Eq. ([^ 

tpn+l = [Mn+lM~^] li)n- (H) 

Or, more generally, using the relation between the wave- 
functions in the m and n layers 

f/’n = tpm, , (12) 


where the exact expression of Tn^m = MnM^^ is known. 
Note that the expression for the transfer matrices is 
disorder-free, which implies that disordered wavefunc¬ 
tion in the first layer propagates into the subsequent lay¬ 
ers just as a clean wavefunction would. This, of course, 
makes sense since only the top layer contains impuri¬ 
ties. A conventional transfer matrix constructed from 
the top surface, however, would always contain impurity 
potentials, and therefore the construction of the whole 
wavefunction would not be as straightforward. 


D. ‘Holographic’ mapping of the impurity potential 


For completeness, we address another question of inter¬ 
est: what is the effective impurity potential experienced 
by an electronic state in the bulk due to the surface impu¬ 
rity. This question can be answered in the same formal¬ 
ism introduced in earlier sections. To compute the effec¬ 
tive impurity potential, we integrate out the first (n — 1) 
layers. The recursion relation for effective potential is 


Vn+i = B 


E-Hn-Vn 


-B^ 


(13) 


with boundary condition Vi = Ump- The ‘holographic’ 
mapping helps us to analytically derive a scaling behavior 
of the effective potential 

Vn=BMn-lM-\ (14) 


where the M^’s are similarly constructed to satisfy the 
Schrodinger equation for the individual layers, Q, and 
to be invertible. However, their boundary condition is 
different from the previous clean case. M„ has to be 
constructed such that the following condition is satisfied: 

Uin,p = BMoM^K (15) 

Since Uimp is a random matrix, it is nontrivial to de¬ 
termine Mq in general. But, we know the object Mn 
propagates just like a clean wavefunction. Therefore, it 
is possible to deduce the scaling of (14,)^ with respect 
to layer index n; which includes the contribution from 
surface impurities as well as clean layers from the top to 
(n — 1)*^ layer. 


III. APPLICATION TO A 2D TOPOLOGICAL 
INSULATOR 

In the previous section, we introduced a general 
transfer-matrix framework for computing the full wave- 
functions of layered systems with surface impurities. In 
this section, a 2D topological insulator modeP^ is em¬ 
ployed to explicitly show how the local density of states 
can be computed in a system with edge impurities. 

Our main results are presented in Fig.[^[^and[^ where 
we use the formalism developed earlier to compute the 
LDOS of the first and second layers of the TI varying the 
disorder strength W. 
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Xj=1 Top layer X|=N 

♦ » # #»♦•#»♦ Layer 1 
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Layer n 
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FIG. 2. A 2D lattice model for a TI with edge impurities. 
Along a:-direction we have a periodic boundary condition that 
a;i = Xn+ 1 , and along y-direction there are N-layers and the 
subscript n indicates the nth layer along the y-direction. Dif¬ 
ferent colors and sizes of dots on the top layer represent the 
random on-site impurity potentials. We study how the ef¬ 
fect of the top edge impurities propagate down into the bulk 
layers. 


A. Model Hamiltonian 

Consider the toy model of a 2D topological insulator. 
In momentum space, 

H{k) = [m — 26(2 — cos — cos ky)]Tz 

+A[txSz sin kx + Ty sin ky] (16) 

where Tt is a Pauli matrix in orbital basis, St a Pauli ma¬ 
trix in spin basis. The lattice spacing is set to a = 1 
such that the momenta kx and ky lie within the inter¬ 
val [—7r,7r]. To introduce an edge state, open boundary 
conditions are introduced in the y-direction, and periodic 
boundary conditions are applied to the a;-direction. The 
intra-layer Hamiltonian and the hopping term between 
layers described in Eq. Q are: 

Hq = [m — 26(2 — cos kx)]Tz + TxSzAsinkx, 

B = hTz- l^Ty. 

The system is in the topological phase if the bands are 
inverted for some range of momentum: sign(mB) > 0. 
For this case, the dispersion of the top and bottom edge 
states are given hy E = zLA sin k^^. 


B. A single layer effective Hamiltonian 

The system is equivalent to a set of parallel ID wires 
coupled by a hopping matrix B in spin and orbital basis. 
We want to construct the matrix which is an essen¬ 
tial building block of a single layer Hamiltonian [Eq. Q] 
and the transfer matrix [Eq. Ill]- To demonstrate the 
method, we construct Mn{kx = 0, E = 0, = 1) here, 

and other {kx,E) cases will be shown in the appendix. 
Note that because the clean Hamiltonian is diagonal in 


momentum and spin space, we only need to analyze the 
orbital space. The Schrodinger equation we need to solve 
is 



,-26 





i’n = 0 , 


where we use ipn+m = P^^ipn- This is one section of 
Schrodinger equation at kx = 0, E = 0, and 5^ = 1. 
There are four transfer eigenvalues and corresponding 
eigenvectors. By taking the determinant of the terms in 
the square bracket, we get pi = Xi, P 2 = 1/Ai, pa = A 2 , 
and p 4 = I/A 2 , where 


Ai 

A 2 


— (m — 26) — y'm? — 4mb + 

A + 2b 

— {m — 26) -I- y/m'^ — Arab + A'^ 

A + 2b 


where Ai ^2 is chosen to be |Ai, 2 | < 1 for 0 < to < 46 
and corresponding eigenvectors are (pi = 4>a = 1+) ^nd 
eigenvectors of 1 /Ai _2 are ^2 = '('4 = |—), where 


1 +) 



I-) 



(18) 


Considering only one spin section, this implies that a 
state |-l-) is localized at the top (n = 1) and a state |—) is 
localized at the bottom (n = N). The interlayer hopping 
operator can be expressed in terms of the eigenvectors: 
B = (6—A/2)|—)(-|-|-|-(6-(-A/2)|-|-)(—|. For this given set 
of eigenvalues and vectors, we can construct an invertible 
Mn in the following manner: 

Mi+isr+i = {p[ - /33)|+)(+| - ip2 - P4)l“)(“li 


where I = n — N — 1. Mn is invertible for n ^ N + 1 and 
satisfies the homogeneous boundary condition Mftf+i = 
0. If TO^ — 4 to6 + <0, eigenvalues are pi = Ae*® 

and Pa = Ae“*® with A = Thus, the effective 

Hamiltonian in the first layer is as follows ([^: 


E - Hf = BMoMp 
= B 


-1 


-N-l 

Pi 


-N-l 

P3 


-N 

Pi 


-N 

P 3 


l+)( + l 


-N-l 

P 2 


-N-l 

Pi 


-N 

P 2 


-N 

Pi 




-)(-| 


(19) 


where N is the number of parallel wires. The effective 
Hamiltonian contains off-diagonal elements only in |±) 
basis. 

Note that the effective single-layer Hamiltonian at 
{kx,E) = (0,0) depends on the width of the system, 
and if a wavefunction contains a component at {kx, E) = 
(0,0), it will also be system-size dependent. Because 
we cannot think of a localized wavefunction dependent 
of the system size for large enough N, we can say no 
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eigenstate localized to an edge sits at {kx^E) = (0,0). 
More relevant Hamiltonian sections at zero energy will 
be kx ^0, which is system size independent in the large 
N limit. More generally, the Hamiltonian sectors not at 
{kx,E) = {27 tI/N, As'mkx) is expressed in the following 
way: 


{+\Hf\+) 

{+\Hf\-) 

{-\Hf\+) 

{-\Hf\-) 


b + Al2 

Xjr-i - 1/ri 

6 + H/2 

l/rs - 1/ri 
b — H/2 p _i -11 

- P3 rs] 


[Pi^/ri- Ps^/rs] , 


1^1 - rs 


[PT^-PI^]- 


( 20 ) 



Here pi^a’s are the eigenvalues of the transfer matrix with 
magnitude smaller than unity, and = {+\(pi)/{—\(j)i)’s 
are the ratios of the overlaps between the transfer-matrix 
eigenstates corresponding to pi with the |-|-) and |—) 
states. We can see that as the {kx,E) approaches to the 
on-shell condition, rp approaches zero and the wavefunc- 
tions have infinitesimal overlap with |—) since (—) 
component is huge. In other words, the Hamiltonian ex¬ 
pressed in this way can be interpreted as a projection to 
the on-shell eigenstates. 

With a set of impurities on the top wire (Fig.[^, to ob¬ 
tain eigenenergies we find energies where the determinant 
of effective single layer Hamiltonian is zero (see Eq. (§). 




FIG. 3. For disorder strength W = 20, histograms of the 
LDOS ratio for different system sizes Nx = 10, 20, 40,80 are 
shown. When system size > 80 the distribution converges 
and the ratio of the LDOS in the first and the second wire be¬ 
comes system size independent, and the average of the LDOS 
ratio in the thermodynamic limit can be estimated. On the 
other hand, system sizes Nx < 40 are in quantum-dot regime 
and not proper to compute the bulk electronic properties be¬ 
cause of their system-size dependence. 


FIG. 4. The LDOS ratio of the first and second layer is 
plotted over the impurity strength W for system size Nx = 
20, 50,120, 200. The edge state at zero energy is primarily 
populated in the first layer from weak to moderate impurity 
strength, W < 5, and then the edge state moves to the sec¬ 
ond and following layers at strong impurity strength. The 
Nx = 200 curve shows the behavior in thermodynamic limit 
as it becomes size-independent, while Nx = 20 curve shows 
quantum-dot behavior for VF > 5, meaning the system size is 
not big enough to see bulk properties. 


In the strong impurity regime we must use a large 
enough system size to correctly see the size independent 
behavior of bulk electronic properties. Here, we distin¬ 
guish the quantum-dot regime from the bulk regime by 
the dependence of physical observables on the system 
size. Figure shows the histogram of the ratio of the 
local density of state in the first and the second wire for 
impurity strength IF = 20 with increasing system sizes 
Nx = 10, 20,40,80. The series of histograms shows size- 
dependence for Nx < 40, the histogram becomes Gaus¬ 
sian shape and size-independent for Nx > 40. Therefore, 
if one wants to numerically obtain physical observables in 
the thermodynamic limit, it is important to use system 
size larger than Nx = 80 for IF = 20 non-magnetic edge 
impurities. The large-size requirement is less stringent 
for smaller impurity strengths, as evident from Fig. 

Figure shows the ratio of the LDOS at zero energy 
with increasing impurity strength. This quantity tells us 
where the edge-state wavefunctions actually reside. A ra¬ 
tio below 1 indicates edge states rooted in the first layer. 
But ratios greater than 1 indicate edge states expelled to 
the the second layer, which are therefore increasingly less 
immune to the disorder. 

The edge states in weak and moderate disorder edge 
state are dominated by the first layer (also see Figure [^. 
This comes hand in hand with a spread of the edge func¬ 
tion Fourier transform: it has broad support away from 
kx — 0 due to impurity scattering. Once the impurity 
strength is comparable to or larger than the band width, 
the edge state is populated less in the impurity layer, and 
























6 


it moves to the second and following layers. All curves in 
Figure]^ show this behavior with a dip at W = 5. While 
we believe that = 200 curve properly describes the 
system-size independent LDOS ratio in thermodynamic 
limit, Nx = 20 curve is only good for VF < 5 and it be¬ 
gins to deviate from = 200 curve for strong impurity 
strength W > 5. 


C. Transfer matrix between single-layer 
wavefunct ions 


Once the wavefunction of the first layer, ipi, is ob¬ 
tained, we next propagate it to the subsequent layers to 
obtain a full profile of the state. This can be done by us¬ 
ing the matrices as in Eq. (12). Let us write down the 
expression of the transfer matrix for the {kx,E) = (0,0) 
case first from layer m to layer n 


= [A”-'"|+)(+| + A" 


|-)(-| 


sm{N -I- 1 — n)9 
sin(A^ + 1 — m)9 ’ 


where = Ae^*® is used as before with A < 1, A' is 
the number of layers. We can see that the |-|-) compo¬ 
nent exponentially decays from the top surface towards 
the other end (n > m), while |—), if it is present, ex¬ 
ponentially increases. Thus, it is apparent that |-|-) is a 
state localized to the top edge, while |—) is localized to 
the bottom edge. However, note that the transfer ma¬ 
trix contains an oscillating term dependent of the system 
size A, just as in the effective single-layer Hamiltonian, 
Eq. (19). It implies that if there is a |-|-) component in 
the wavefunction at exactly {kx,E) = (0,0), its oscil¬ 
lating part is dependent on the number of layers, which 
doesn’t make sense in the physical picture where A is 
much larger than the localization length of edge state. 
Therefore, we can say |-|-) component at {kx,E) = (0,0) 
must be vanishingly small as the system size is increased. 

The transfer matrix from the first layer to the layer 
for a general (kx,E) is expressed in the following way: 


(+|T„^i|-|-) — 


1 


[prVi-prVa] 


ri - rs 

ri - rs 

= i/r. -T/n ■ '’’■'''’'=1 ■ 


( 21 ) 


with Vi = {+\cl),)/. 

We apply this transfer matrix to the first layer wave- 
function to obtain wavefunctions in the bulk layers. In 
Fig.m the disorder-averaged probability density in mo¬ 
mentum and layer basis P{kx,n) = \'4’n{kx)\'^ is plot¬ 
ted for shown impurity strength W. In the weak disor¬ 
der regime, W = 0.1, where impurity strength is much 
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Layer n 
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Layer n 


2 4 6 8 10 12 14 

Layer n 



FIG. 5. The disorder-averaged LDOS at A = 0 in wire index 
n and momentum space kx of edge state at zero energy is 
plotted for increasing impurity strength. Starting from the 
clean edge state where the LDOS is only at kx = 0, the LDOS 
is spread out in momentum space then shifted to the second 
layer. The insets are showing the fc 3 ;-integrated LDOS as a 
function of layer n. Large enough system size Nx = 180 is 
chosen such that the averaged LDOS is size-independent. 500 
disorder realizations are averaged. 


smaller than the energy gap, the probability density is 
concentrated near kx = t) at zero energy. As the impu¬ 
rity strength increases the probability density gains width 
in momentum space and its weight is shifted to the sec¬ 
ond layer. While this trend is quite strong already with 
W = 10, in strong impurity regime W = 40—which is 
much larger than the bandwidth—the zero energy wave- 
function is completely absent in the first layer, but occu¬ 
pies the subsequent layers in a narrow range of momen¬ 
tum space. This indicates that the wavefunction has been 
pushed to the next layer and behaves as if the system is 
clean. 

This behavior of the local density of states is shown for 
non-magnetic edge impurities, which cannot affect the 
transport properties of helical edge states in 2D topo¬ 
logical insulators. Therefore, the modification of LDOS 
should be discussed separately from the change of trans¬ 
port nature, at least in 2D. For a strip geometry, the 
transport is studied for both non-magnetic and magnetic 
edge impurities in the following sections. 

The widths of probability density P{kx, n) = |'!/'n(fca;)P 
in the first and second layer at zero energy are plotted 
in Fig. The width in momentum space is indicative of 
how disordered the edge state is due to impurities. In the 
weak and strong disorder limit the wavefunction behaves 
like a clean system in the LDOS shape as shown in Fig. 
With strong disorder, the width of P{kx, n = 1) saturates 
near 0.3, although it carries little weight in that limit: 
P{kx, n = 1) ^ 1. Meanwhile, the width of P{kx, n = 2) 
increases and decreases again as the impurity strength is 
varied, peaking at around the bandwidth of the system 
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FIG. 6. For system size = 120, the width of LDOS in the 
first and the second wire in momentum space kx is plotted for 
as a function of disorder strength W for non-magnetic impu¬ 
rities. In clean system limit, the edge state LDOS P{kx) is 
concentrated at fca, = 0 and its width is infinitesimally small. 
As the disorder strength increases, the LDOS spreads in mo¬ 
mentum space. However, for the stronger disorder, the first 
layer LDOS takes on all possible momenta and its width sat¬ 
urates, while the second layer LDOS become concentrated 
around kx = 0 again. 


(iy-8). 


D. Magnetic edge impurities 

The same calculation was repeated for a system with 
magnetic edge impurities. We simply needed to extend 
the Hamiltonian to have two spin-sections and intro¬ 
duce random magnetic impurities, V{xi) = Vi • s, where 
three component random variable Vi = {Vf, , Vi). We 
found that to simulate the bulk regime for IT = 20 the 
system size needs to be at least = 400 as opposed to 
Nx = 120 for the non-magnetic edge impurity case. In 
other words, the lower bound of the system size to see the 
thermodynamic properties is much larger and it becomes 
computationally challenging even for 2D system. Fig- 
urej^shows the width of the LDOS distribution P{kj;, n) 
for n = 1,2. The data is qualitatively similar to the 
non-magnetic case, which demonstrates the universality 
of the result. 


IV. TRANSPORT BEHAVIOR 


FIG. 7. LDOS momentum-space width in the first and the 
second layers with magnetic impurities, at system size Nx = 
120. The overall behavior is the same with the case of non¬ 
magnetic edge impurity case (Fig.[^, but much larger systems 
are required to see the size-independent physical properties in 
bulk regime. 


non-magnetic impurities, edge modes can not backscat- 
ter, and their conductance remains quantized at the value 
of the clean system, despite the local density of states 
associated with it changing its support between the lay¬ 
ers. To clarify the transport nature of the system with 
magnetic and non-magnetic impurities along the edge, in 
this section we study conductance of the systems using 
Landauer-Biittiker method. 

Imagine a system where two semi-infinite leads are con¬ 
nected to a disordered region at the ends x = 1 and 
X = Nx- Landauer and Bii” ttikei!22H2i! related the con¬ 
ductance with the transmission coefficient through the 

2 

disordered region: g = \Tni, where g is conductance, 
Tjvi is a transmission coefficient through the disordered 
region from site 1 to site N. Using linear response for¬ 
malism, Fisher and Le^^ expressed the transmission co¬ 
efficient in terms of Green’s functions: 



^ lGni^ rG 


t 

N1 


( 22 ) 


where F^ = Gni is a Green’s function from 

site 1 to iV renormalized by the presence of the leads, 
and El is a self-energy of the semi-infinite left lead. Each 
term in this formula can be computed recursively, such 
that the conductance of a long system can be obtained 
with a reasonable computation effort. A good review of 
the detailed calculation can be found in Ref. [26l 


The local density of states discussed in the last sec¬ 
tion can be probed by angle-resolved photoemission 
spectroscopy.® The transport along the edge states, how¬ 
ever, is provides additional information independent of 
the local density of states. For instance, in the case of 


A. Non-magnetic impurities case 

Gonsider the 2D topological insulator system intro¬ 
duced earlier with non-magnetic impurities along the top 




















edge. Because the Hamiltonian is diagonal in spin basis 
without magnetic impurities, we can consider transport 
in just one spin sector. When the chemical potential is in 
the energy gap, neither backscattering nor scattering into 
the bulk is possible. Therefore, the conductance must re¬ 
main quantized even in the presence of edge impurities. 
This is indeed what our calculation shows. Indeed, when 
the disorder is non magnetic, the transport behavior does 
not reflect the development of the LDOS. 


B. Magnetic impurities case 

The two opposite-spin, counter-propagating, chiral 
edge modes couple as soon as magnetic edge impurities 
are introduced. As a result, transport through the dis¬ 
ordered edge is suppressed, while the transport through 
the clean edge remains unaffected. Therefore, we expect 
the total conductance to rapidly approaches e^/h when 
introducing and increasing magnetic edge disorder. 

However, in the strong-impurity quantum-dot limit, 
W S> we found that the conductance recovers its 
clean system value of 2e^//i. In this regime the impuri¬ 
ties are strongly bound to electrons at energies far away 
from the Fermi energy, and they play negligible roles in 
the transport at the Fermi energy. Put another way, 
strong disorder pushes the edge modes to the next layer 
where they effectively become weak scatters. This be¬ 
havior is clear in Fig.[^ which shows the conductance vs. 
the disorder strength for different system sizes. We can 
see that in the intermediate range of impurity strength, 
the conduction through the disordered top edge is sig- 
nihcantly suppressed due to magnetic impurities, while 
the conductance recovers up to 2e^//i value in the strong 
impurity limit. We note that in the thermodynamic limit 
Nx —>■ 00 , the conductance is always jh for any disor¬ 
der strength, illustrated in Fig.[^ Our calculation reflects 
the non-monotonic dependence of the localization length 
of the scattered edge mode on disorder strength, which 
is consistent with the expulsion of the LDOS from the 
disordered first layer. 


V. SUMMARY AND CONCLUSION 

The effects of surface disorder on systems with sur¬ 
face states are at the focus of our work. The interest in 
this problem rose after an exact-diagonalization analysis 
showed that the response of a 3D topological insulator 
to surface disorder is non-monotonic: First the surface 
states become diffusive, and their conductance is sup¬ 
pressed. At a finite disorder strength, however, the sur¬ 
face states mean free path recovers, and they reconstitute 
at the disorder-free second layer of the TP^. 

In this manuscript we developed a formalism that re¬ 
duces solving a bulk D-dimensional Hamiltonian to a 
surface-only diagonalization problem. This, in principle, 
enables an exact-diagonalization analysis of much bigger 



FIG. 8. Edge channel conductance vs. impurity strength 
when magnetic edge impurities are present. Three system 
sizes Nx = 20, 50,100, 200 are used. The edge channel con¬ 
ductance drops to e^/h, initially. The conductance recovers 
to 2e^//i as the system enters quantum-dot regime at strong 
magnetic impurity strength. This reflects the nonmonotonic 
behavior of the localization length of the edge mode as disor¬ 
der increases. 



FIG. 9. Disorder averaged conductance with increasing sys¬ 
tem size Nx. for three impurity strength W. The conductance 
exponentially decreases with system size as anticipated in the 
magnetic edge impurity case. 


bulk systems than previously possible. Roughly speak¬ 
ing, our method constitutes a systematic integrating out 
and elimination of the bulk degrees of freedom layer by 
layer. We show how to carry this procedure out with a 
technique reminiscent of transfer-matrix methods. From 
the resulting surface-only diagonalization problem we are 
able to reconstruct the wavefunctions of the bulk layers 
and study the effect of disorder on them. Our method is 
generally valid for any system composed of coupled lay¬ 
ers but becomes particularly useful for topological sys¬ 
tems whose surface states are fundamentally dependent 
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on the existence of the gapped bulk. 

We used our method to study 2D topological insulators 
with edge impurities. As the strength of edge impurities 
increases, we found that edge states near zero energy 
become more concentrated up to the moderate impu¬ 
rity strength (i.e., comparable to the bandwidth), and 
then the edge states are gradually pushed to the second 
layer. The ratio of the local density of states in the first 
two layers in Figure shows this behavior for different 
system sizes. Furthermore, the width of the edge state 
in momentum space in the second layer reaches a max¬ 
imum at moderate impurity strength and then narrows 
at stronger impurity strength. Despite the strong disor¬ 
der, the edge states have momentum restored to being 
a good quantum number. This non-monotonic response 
to the surface disorder is equally true for magnetic and 
non-magnetic disorder (see Figures and [^. 
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The transport properties, however, show a sharp con¬ 
trast between the magnetic and non-magnetic cases: non¬ 
magnetic edge impurities does not affect the edge trans¬ 
port properties, while magnetic edge impurities immedi¬ 
ately induce a finite localization length due to backscat- 
tering, and initially suppress the edge conductance (see 
Figure |^. In our numerical simulations that include 
magnetic impurities, however, we found that the con¬ 
ductance then recovers to the clean-system values when 
the impurity strength exceeds the system size. 

This finding does not imply that the edge states com¬ 
pletely decouple at large but finite disorder. Rather, it 
is an indication that the localization length of the spin- 
orbit locked modes exhibits a non-monotonic localization 
behavior as a function of disorder. In our simulations, 
perfect conductance will be recovered when the localiza¬ 
tion length exceeds the system width. This is precisely 
the regime we nicknamed the quantum-dot regime. The 
localization length inferred from transport calculations 
must be proportional to the inverse of the average width 
of the edge LDOS in momentum space, for momenta par¬ 
allel to the edge. Therefore, our transport results are yet 
another manifestation of the reconstitution of the edge 
mode at the second layer, which is disorder free. In the 
limit of infinitely strong disorder, we expect the localiza¬ 
tion length of the reconstituted edge state to diverge and 
translational invariance is regained. 


In future work, we intend to apply our method to the 
3D topological insulator with surface disorder. As we em¬ 
phasize throughout, one must use large enough systems 
in order to explore the bulk properties in thermodynamic 
limit (see Figure|^and|^. This makes 3D topological sys¬ 
tems with moderate to strong surface impurities harder 
to study using the ED as the cost of calculation increases. 
Our analytic approach will then be very useful, since it 
reduces the computational cost dramatically. 
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APPENDIX: DERIVATION OF THE EFFECTIVE where j = 1,2 and / = 3,4. Then, we can think the 

HAMILTONIAN AND THE TRANSFER MATRIX following construction of the object M: 


The construction of the object M in Sec. |TT] is a cen¬ 
tral element for further computation in any examples. In 
this appendix, we show how the object M in Sec. |III| was 
constructed for given momentum kx and energy E. We 
work in the momentum space since the clean Hamilto¬ 
nian is diagonalized in momentum space. But, as the 
layer degree of freedom is integrated out, the effective 
Hamiltonian is dependent of energy nontrivially as can 
be seen in the recursive Eq. 0. Therefore, the disper¬ 
sion relation cannot be immediately deduced from the 
single layer effective Hamiltonian. Rather, we compute 
the quantity det[E — — Pimp] in momentum space as 

varying energy E such that the determinant is zero. In 
this way we find eigenenergies and then eigenstates of a 
surface disordered system. 

The object M satisfies the layered Schrodinger equa¬ 
tion Eq. Q. Thus, the construction is convenient in 
terms of the wavefunctions of the Schrodinger equation, 
Eq. 0. The most general expression will be 

Mi+N+i = ^ (23) 

i,i=l,2,3,4 

with the vanishing boundary condition at the last layer 
Mtv+i = 0. There are different ways to build M and we 
introduce one way for on-shell E = A sin kx and off-shell 
E ^ A sin kx states. 


1. On-shell states, E = Asink^ 


Considering only one spin sector of the Hamiltonian 
as in Sec. |HIB[ the eigenvalues and eigenvectors of the 
transfer matrix over the layer are straightforwardly ob¬ 
tained from the Schrodinger equation with the replace¬ 
ment Ipn-l = -pi’n and 'tpn+l = pi>n, 

B^ij;nE[Ho{kx) - E]’4}n +pipn = (24) 

Solving a 2x2 matrix equation, we obtain four eigenvalues 
Pi=i, 2 , 3,4 with eigenvectors (/i'i=i, 2 , 3 , 4 - For on-shell states, 
we immediately know that two eigenvectors are |-|-), 


<('1 



1 

71 




with eigenvalues \pi=^^ 4 \ < 1 so that <(>j= 3,4 are physi¬ 
cally localized states at the top layer n = 1. This is to 
satisfy the vanishing boundary condition of the clean sys¬ 
tem with two edge-localized wavefunctions. To construct 
the object M, we need not only the vanishing boundary 
condition, but also M must be invertible. To do that, 
consider the decomposition of 4 >i=i ^2 into |-|-) and |—): 


Pj = l,2 — 


■3 1 _ ^3 


v^(l 


72 


1 


(25) 


From this set of eigenvectors, we can construct two copies 
of |—)’s behaving differently over the layers. Specifically, 

1 A 

|-(h0) = -^^3 - 


Mn+1 — [l“(l,3)) - |-(2,4))] (-| + [<('3 - 4’i] ( + 1 

which is zero. For n G [0,iV], is of course non-zero 
as eigenvalues of the eigenvectors are different, and 
is generally invertible. For nth layer. 


(+|M„|+)= pn-N-l_pn-N-l^ (27a) 

{+\Mu\-)= 

I ''^2 r n—N—1 11 /'07K1 

+ ^ P4 -P2 (27b) 

n>2 

(-|M„|+) = 0 (27c) 

(-|M„|-)= - pn-N-i (27d) 


From this, the construction of the transfer matrix and the 
effective single layer Hamiltonian is following: T„<_m = 


(+|7n.,_m|+) 


{+\Tn<-m\ — } 


n—N—1 n—N—1 

P 3 _P 4 _ 

m—N—1 m — N—1 

P 3 P 4 

n—N—1 _ A 2 „n—N — l 
Bi Pi _B 2 _ 

^m—N—l nn—N—1 

Pi -P 2 


(28a) 


^ m-Af-l _ Ai ^m-N-1 

\P3 P 4 ) \ B2 P'2 Bi Pi 


(P3 -Pi )(Pl -P2 ) 

(28b) 


(~|7n<-m|+) — 0 




n—N—1 n—N—1 

Pi _P2_ 

m—N—1 m — N—1 

Pi P 2 


(28c) 

(28d) 


Expressing B = {b + A/2)|-|-)(-| -I- (6 - A/2)|-)(-|-|, the 
effective single layer Hamiltonian is: 


{+\E-Hf\+) = 

(6-A/2)(-|ro, 

-il+) 

(29a) 

{+\E-Hf\-) = 

(6-A/2)(-|ro. 

-il-) 

(29b) 

{-\E-Hf\+) = 

{b -l- A/2)(-|-|ro<- 

-i|+) 

(29c) 

{-\E-Hf\-) = 

(6 -l- A/2)(-|-|ro<- 

-i|+) 

(29d) 


One can see that the effective Hamiltonian is still system- 
size dependent as Hf^{kx = 0,E = 0) case computed 
in the main manuscript. Plus, E — E[f^ contains no 
zero eigenvalues, implying that any finite system cannot 
have the energy dispersion E = A sin kx due to finite- 
size effects. Therefore, the expression of the Hamiltonian 
for on-shell states is not useful for actual computation. 
Rather, we need the Hamiltonian expression of off-shell 
states for finites size system with surface impurities, for 
which we find an analytic expression of the effective single 
layer Hamiltonian and let A —>■ oo for the semi-infinite 
limit. 


(26) 
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2. OfF-shell states, E 7 ^ Asinfej, 

For a system with surface impurities, eigenstates are not described by the clean system dispersion E = Asink^^ 
rather, the state has a mix of different momentum components in each given energy. More concretely, if the size of the 
syst em along the periodic boundary condition is Nx, the exact edge state dispersion of the clean system discussed in 
Sec. |III| is E = As\a{2TTl/Nx) with integer 1. Only at those discrete set of energies, we have two eigenvectors parallel 
to 1+) and the discussion in the previous section applies. Except those on-shell points, we have the following general 
set of eigenvectors. 


(j)i = 






(30) 


with eigenvalues ^ 7 = 1 , 2 , 3 , 4 - Without loss of generality, let us say |pj= 2 , 4 | > 1- Each eigenvectors can be written as 
the sum of |-|-) and |—) like Eq. (251. Then, we similarly construct two pairs of |-|-) and |—) by the superposition of 

^j = l,2,3,4- 


l+(j.o) - 


1 


1 


Al _ A \B 


1 


I (J.O) Bj Bi [ 


1 

A 


Next, we can construct the object satisfying the vanishing boundary condition at n = N" + 1. 


(31) 


Mn+I — (|-(1,2)) - |-(3,4)))(-| + (l + (l,2)) - l + (3,4)))( + l (32) 

where we intentionally split (/>2 and ^4 in each term so that no terms vanish in —>■ 00 limit. Explicitly, the 

components are 


1 


Ai j A 2 j 




1 


^3 _ 

B3 B4 


(-|-|M/+Ar+i|-|-) — ^ ^ 

Bi B 2 

(-|M/+Ar+l|-k) = ^ ^ [p[ - P 2 ] - A 3 _ Ai [P3 - pi] 

Bi B 2 S 3 S 4 

Bi B2 [^1 ~ ^2] ~ B3 Bj [P 3 ~ Pi] 

Ai A2 As A4 

1 


^ J _ 

b/\ 


{+\Mi+n+i\-) — 



B2 1 

1 


Bi , 

[ a /1 


S3 _ B^ 

As A4 


a/^\ 


(-|M/+Ar+i|-) = ^ ^ 

Ai A 2 

We are interested in the behavior of edge states near the top (n = 1). Thus, in the limit of N" —)■ oo, 

and the terms with ( ^ dominates. 

\Pj = 2,i J 


(33a) 

(33b) 

(33c) 

(33d) 

= n—A ^—1 —)■ —oo 


(+|T'„.i-m|+) — 


(~|Tri<-m|+) — 


A^ n—m A 3 ^n—m 

Si /^1 S 3 ^3 

Ai _ As ’ 

Si S3 

n—m 71—m 

Ps ~ Pi 

Bi _ B 3 ’ 

Ai A3 


(+|2n.i-m|~) — 




n—m n—m 

Pi Pz 

Ai _ As ’ 

Si S3 

S3 n—m _ Si n—m 

A^Ps A^Pl 

Si _ ^ 

Ai As 


(34) 


And the expression for the Hamiltonian is following Eq. (29a|-(29d). 
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